#################
#Estimation procedures for spatial rebels
#Project: Spatial Rebels
#Authors: NMJW
#Data created: 17 May 2013
#Last Update: 17 April 2014
#################
#  SYSTEM REQUIREMENTS
#  This script was developed and run under R64.app operating under WIN7.
#  R version: R 2.15.0
#################


###############
# Naive Log-Likelihoods
###############
n.weibull <- function(theta, k, X, Y){
	beta <- theta[1:k]
	esigma <- exp(theta[k+1])
	u <- (1/esigma)*(Y - X %*% beta)
	logl <- u - exp(u) -log(esigma)
	return(sum(logl))
}
n.loglogistic <- function(theta, k, X, Y){
	beta <- theta[1:k]
	esigma <- exp(theta[k+1])
	u <- (1/esigma)*(Y - X %*% beta)
	logl <- u +log((1 + exp(u))^(-2)) -log(esigma)
	return(sum(logl))
}
n.lognormal <- function(theta, k, X, Y){
	beta <- theta[1:k]
	esigma <- exp(theta[k+1])
	u <- (1/esigma)*(Y - X %*% beta)
	logl<- -.5*log(2*pi) - .5*log(esigma^2) -.5*(u*u)
	return(sum(logl))
}


###############
# Spatial Log-Likelihood functions
###############
sp.weibull <- function(theta, k, X, Y, Wy, eig){
	beta <- theta[1:k]
	rho <- theta[k+1]
	if (rho < -1 | rho > 1) return(NA)
	esigma <- exp(theta[k+2])
	u <- (1/esigma)*(Y - rho*Wy - X %*% beta)
	logl <- log(1 - rho*eig) -log(esigma) + u - exp(u)
	return(sum(logl))
}
sp.loglogistic <- function(theta, k, X, Y, Wy, eig){
	beta <- theta[1:k]
	rho <- theta[k+1]
	if (rho < -1 | rho > 1) return(NA)
	esigma <- exp(theta[k+2])
	u <- (1/esigma)*(Y - rho*Wy - X %*% beta)
	logl <- log(1 - rho*eig) -log(esigma) +u +log((1 + exp(u))^(-2))
	return(sum(logl))
}	
sp.lognormal <- function(theta, k, X, Y, Wy, eig){
	beta <- theta[1:k]
	rho <- theta[k+1]
	if (rho < -1 | rho > 1) return(NA)
	esigma<- exp(theta[k+2])
	u <- (1/esigma)*(Y - rho*Wy -  X %*% beta)
	logl <- log(1 - rho*eig) -.5*log(2*pi) - .5*log(esigma^2) -.5*(u*u)
	return(sum(logl))
} 


###############
# Estimator Functions
###############


# naive duration
nduration <- function(fmla, dat,dist){
	# prepare data from formula
  fmla <- as.Formula(fmla)
  dat <- model.frame(fmla,data=dat)
	Y <- log(as.matrix(model.part(fmla, data=dat, lhs = 1)))
	Y.org <- as.matrix(model.part(fmla, data=dat, lhs = 1))
	N <- length(Y)
	X <- model.matrix(fmla, data = dat, rhs = 1)
	k <- ncol(X)

	# Define initial values for the parameters
  theta.start <- c(rep(0, k), 1)
  names(theta.start) <- c(colnames(X), "log(scale)")

  # Calculate the maximum likelihood
  mle <- maxLik(dist, method="BFGS", start=theta.start, k=k, Y=Y, X=X)
  par <- mle$estimate
 	
  # expected durations
  yhat <- NULL
  lnyhat <- NULL
	sigma <- exp(par[k+1]) 
 	mu <- X %*% par[1:k]
 	# weibull
 	if (deparse(substitute(dist))=='n.weibull') {	
 		yhat <-  exp(mu) * gamma(1+ sigma)
 		lnyhat <- mu -digamma(1)
 	}
  # loglogistic
 	if (deparse(substitute(dist))=='n.loglogistic') {	
 	  if (sigma<1) {   
      yhat <- exp(mu)*(pi*sigma) / sin(pi*sigma)    # loglogistic
 	  } else {
      yhat <- NA
      print('Warning: expected duration not defined for scale < 1')
 	  }
 	  lnyhat <- mu
 	}
	# lognormal
	if (deparse(substitute(dist))=='n.lognormal') {	
		yhat <- exp(mu + .5*sigma^2)
		lnyhat <- mu
	}
	
  out <- list(mle=mle, yhat=yhat, lnyhat=lnyhat,N=length(yhat),y=as.matrix(model.part(fmla, data=dat, lhs = 1)),data=dat,X=X,Y=Y,par=par,Y.org=Y.org)
  return(out)
}


###############
###############
predict.duration <- function(dist,X=X,par=par){
	 # expected durations
k <- ncol(X)
  yhat <- NULL
  lnyhat <- NULL
	sigma <- exp(par[k+1]) 
 	mu <- X %*% par[1:k]
 	# weibull
 	if (deparse(substitute(dist))=='n.weibull') {	
 		yhat <-  exp(mu) * gamma(1+ sigma)
 		lnyhat <- mu -digamma(1)
 	}
  # loglogistic
 	if (deparse(substitute(dist))=='n.loglogistic') {	
 	  if (sigma<1) {   
      yhat <- exp(mu)*(pi*sigma) / sin(pi*sigma)    # loglogistic
 	  } else {
      yhat <- NA
      print('Warning: expected duration not defined for scale < 1')
 	  }
 	  lnyhat <- mu
 	}
	# lognormal
	if (deparse(substitute(dist))=='n.lognormal') {	
		yhat <- exp(mu + .5*sigma^2)
		lnyhat <- mu
	}
	
  out <- list(yhat=yhat, lnyhat=lnyhat,N=length(yhat))
  return(out)
	
	}


# statial duration
#spduration <- function(fmla, dat, W, dist){
spduration <- function(fmla,uniqID,dat, W, dist){
	#HERE A TEST IF UNIQUE ID IS REALLY UNIQUE	
	# prepare data from formula
  fmla <- as.Formula(fmla)
  fmlab <- as.character(fmla)
  fmla.frame <- as.Formula(paste(fmlab[2],fmlab[1],fmlab[3],"+",uniqID,sep=""))
  dat <- model.frame(fmla.frame,data=dat)
  ##NEED TO IMPOSE ORDER ON DAT AND W
   want <- which(colnames(W) %in% dat[,dim(dat)[2]]) 
   W <- W[want,want]
	X <- model.matrix(fmla, data = dat, rhs = 1)
	Ys <- as.matrix(model.part(fmla, data=dat, lhs = 1))
	Y1 <- log(as.matrix(Ys[,1]))
	Y2 <- log(as.matrix(Ys[,2]))
	Y.org <- Ys[,2]
	Y <- Y2	
	N <- length(Y1)
	k <- ncol(X)
 
  
	# calculate spatial lag based on lower half of Y vector
	YY <-  rbind(Y1,Y2)
	WY <- as.matrix(W%*%YY)
	Wy <- WY[(N+1):(2*N)]
	
	# identity matrix
	I <- matrix(nrow=N, ncol=N, 0)
	diag(I) <- 1
	
	# W is 2N x 2N; only lower right relevant for determinant
	W <- W[(N+1):(2*N),(N+1):(2*N)]
	
	# eigenvalues
	eig <- eigen(W, only.values=TRUE)$values	
    
	# Define initial values for the parameters
  theta.start <- c(rep(0, k), 0, .5)
  names(theta.start) <- c(colnames(X), "rho", "log(scale)")

	# Calculate the maximum likelihood
	#mle <- optim(theta.start, dist, k=k, Y=Y, X=X, Wy=Wy, eig=eig, hessian=T, method="BFGS")
	mle <- maxLik(dist, method="BFGS", start=theta.start, k=k, Y=Y, X=X, Wy=Wy, eig=eig)
	par <- mle$estimate
	
	# expected durations
	yhat <- NULL
  lnyhat <- NULL
  mu <- solve(I-par[k+1]*W) %*% (X %*% par[1:k])
  sigma <- exp(par[k+2]) 
  # weibull
 	if (deparse(substitute(dist))=='sp.weibull') {	
      yhat <-  exp(mu) * gamma(1+ sigma)
      lnyhat <- mu -digamma(1)
 	}
	# loglogistic
	if (deparse(substitute(dist))=='sp.loglogistic') {	
	  if (sigma<1){   
	    yhat <- exp(mu)*(pi*sigma) / sin(pi*sigma)
	  } else {
	    yhat <- NA
	    print('Warning: expected duration not defined for scale < 1')
	  }
	  lnyhat <- mu
	}	
	# lognormal
	if (deparse(substitute(dist))=='sp.lognormal') {	
	    yhat <- exp(mu + .5*sigma^2)
      lnyhat <- mu
	} 

  out <- list(mle=mle, yhat=yhat, lnyhat=lnyhat,N=length(yhat),data = dat,X=X,Y=Y,I=I,par=par,W=W,Y.org=Y.org)
  return(out)
}




###############
###############
predict.sp.duration <- function(dist,X=X,I=I,par=par,W=W){
k <-ncol(X)
yhat <- NULL
  lnyhat <- NULL
  mu <- solve(I-par[k+1]*W) %*% (X %*% par[1:k])
  sigma <- exp(par[k+2]) 
  # weibull
 	if (deparse(substitute(dist))=='sp.weibull') {	
      yhat <-  exp(mu) * gamma(1+ sigma)
      lnyhat <- mu -digamma(1)
 	}
	# loglogistic
	if (deparse(substitute(dist))=='sp.loglogistic') {	
	  if (sigma<1){   
	    yhat <- exp(mu)*(pi*sigma) / sin(pi*sigma)
	  } else {
	    yhat <- NA
	    print('Warning: expected duration not defined for scale < 1')
	  }
	  lnyhat <- mu
	}	
	# lognormal
	if (deparse(substitute(dist))=='sp.lognormal') {	
	    yhat <- exp(mu + .5*sigma^2)
      lnyhat <- mu
	} 

  out <- list(yhat=yhat, lnyhat=lnyhat,N=length(yhat))
  return(out)		
}

